Minimum Energy Configurations of Repelling Particles in Two Dimensions 



iJ 



OO 



E. A. Jagls 
Centra Atomico Bariloche 
Comision Nacional de Energia Atomica 
(8400) S. C. de Bariloche, Rw Negro, Argentina 

Geometrical arrangements of minimum energy of a system of identical repelling particles in two 
dimensions are studied for different forms of the interaction potential. Stability conditions for the 
triangular structure are derived, and some potentials not satisfying them are discussed. It is shown 
that in addition to the triangular lattice, other structures may appear (some of them with non-trivial 
unit cells, and non-equivalent positions of the particles) even for simple choices of the interaction. 
The same qualitative behavior is expected in three dimensions. 
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I. INTRODUCTION 

Nature teaches us that crystalline structures —namely, 
the periodic spatial arrangement of atoms— are the min- 
imum energy configurations (MECs) of systems with a 
great number of particles (at least in the case where these 
particles come in only a small number of different types). 
The determination from first principles of the most stable 
crystalline structure of a substance, given the properties 
of its constituent atoms and their interactions is a com- 
plicated minimization problem for which exact methods 
do not exist. A usual way of determining the MECs is 
to compare the energy of different proposed structures, 
and pick up the lowest energy one. The numerical simu- 
lation of finite systems, using the technique of simulated 
annealing may be of great help in this process, since if 
the cooling down of the system is sufficiently slow during 
the simulation, then we expect that the particles accom- 
modate to its MEC. 

When all particles in the system are identical and may 
be considered as point-like, interacting through a poten- 
tial energy that depends only on their relative positions, 
the problem simplifies greatly. In real pure materials, 
most of the MECs correspond to the high symmetry 
structures hep, bcc, fee, and sc. Other structures (no- 
toriously the phases of carbon and ice) are largely due to 
the directionality of atomic orbitals. We will concentrate 
here in the case of isotropic interactions. It is known 
that even in this rather complicated radial inter- 

action potential U (r) (for instance, an oscillating poten- 
tial) may give, rise to complex structures (quasicrystals, 
for instance). El If we restrict to the case of repelling in- 
teractions {U' (r) < 0), with only an external pressure 
preventing the particles to move away from each other, 
the MECs seem to be limited to the above mentioned 
hep, bcc, fee, and sc structures. In the case of particles 
in two dimensions, and under the conditions that all par- 
ticles are identical and the interactions are repulsive, the 
general believe is that the triangular structure (TS) is al- 
ways the MEC (I will use 'triangular' since it is the most 
commonly used term, although from symmetry consider- 



ations the word 'hexagonal' would be more appropriate). 

The aim of this work is to analyze the possible exis- 
tence of structures other than the triangular, for identi- 
cal particles interacting in two dimensions. The case of 
power-law, cut-off power-law, and potentials with a hard- 
core plus a soft repulsive shoulder are discussed in detail. 
It will be shown that the TS may not be the MEC even 
for some 'simple' forms of U (r). 

The paper is organized as follows. In the next sec- 
tion the local stability of the TS against a displacement 
of a single particle is discussed. The MECs of a fam- 
ily of potentials that does not satisfy this criterium are 
shown in Section III. In Section IV, a stability criterium 
against a global deformation is derived, and some po- 
tentials that do not satisfy it are identified. Finally, in 
Section IV there is a short summary and some discussion 
about topics that are not deeply considered in the paper, 
namely the case of three dimensional systems, the effect 
of temperature and other dynamical effects. 



II. LOCAL STABILITY OF THE TRIANGULAR 
STRUCTURE 

We will consider a system of identical particles lying 
on the plane, and interacting through a two-body repul- 
sive spherical potential U (r), U' (r) < 0, where r is the 
separation between particles. The system is supposed to 
be constrained by an external pressure P. We will try 
to find the geometrical configuration of the particles that 
minimizes the free energy of the system. The MEC must 
be obtained by minimizing the enthalpy H = E + PV, 
where V is the volume, and the energy E is given by 
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C/(|r,;-r,|). 



(1) 



The complete solution of this problem by analytical 
means is out of our possibilities. But since the TS is 
our initial guess to the MEC, we can start by analyzing 
its local stability. 
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Clearly, if the TS is the MEC, the energy of the system 
must increase when a single particle is slightly displaced 
from its lattice position. The potential around a given 
site (taken to be the origin of coordinates) created by a 
particle at a generic position Tq is, up to second order, 

of the form U" (ro) Sx^ + ^^^6y^ + U' (ro) dx + U{ro), 
where dx (Sy) is the coordinate of the tested point mea- 
sured along (perpendicular to) Tq. This potential must 
be summed up for all particles, and for lattices with ro- 
tational symmetry C3 or higher (as it is the case for the 
TS) it must reduce to a isotropic form. Considering the 
invariance of the trace of cuadratic forms under rotations, 
the cuadratic part of the final effective potential can be 
written as U2 (Sr/a) /2, with 
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(2) 



where the sum is over all particles of the lattice. The in- 
troduced parameter a is arbitrary, but it will be taken to 
be the lattice parameter of the TS, in such a way that U2 
can be directly compared for lattices with different lat- 
tice parameters. The positiveness of U2 is the condition 
for the stability of |t|he lattice under small displacements 
of a single particle.! 
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FIG. 1. The functions f (7) (thick lines) and f (7) (thin 
lines) to be used for the calculation of the rigidity of the trian- 
gular lattice against the displacement of a single particle and 
a global rescaling (see text). Dashed lines are the contribution 
of nearest neighbors only. 

Let us consider some particular cases. For U (r) = 
Uq {ro/ry all terms of the sum in (^ are positive and 
the sum is convergent for all 7 > 0. The value of U2 is 
given by 

U2 = fi (7) Uo {ro/ay , (3) 
with the adimcnsional function /i (7) given by 

/i(7)=7'E(^''/«)"'"'- (4) 



The function /i (7) is plotted in Fig. |^. In this fig- 
ure, the contribution to /i coming only from nearest 
neighbors is also shown, and it is clear that the con- 
tribution of particles at larger distances becomes more 
relevant as 7 goes to zero. We note in addition that 
U2 vanishes for 7 — > 0. This is related to the fact that 
7 — > gives a logarithmic interaction between particles 
(lim..y^o+ ^ = ~lii(a;)), that corresponds to two- 
dimensional charges, and it is well known that the equi- 
librium configuration in this case has all charges at the 
borders of the system, so no TS is stable in this case. 

Now we will turn to cases when U2 can be negative. 
First of all we should notice that (^ is proportional to 
the Laplacian of U (r), namely 



U2 



^E A{/(r,). 



(5) 



so the problem has an electrostatic analog. Considering 
the electrostatic problem AU = — p, if we look for poten- 
tials U (r) such that AU < for all r, we need a positive 
(or zero) charge density p at all distances. But if in ad- 
dition we require that U (r) vanishes sufRciently rapid 
for r — > 00, we are forced to locate a negative charge 
at the origin (of the same absolute value than the inte- 
grated positive charge). This choice produces a poten- 
tial U (r) that goes to — cxd at the origin, i.e., it would 
not be repulsive at short distances. This shows that U2 
cannot be negative at all distances for repulsive short 
range interactions. However, since expression must 
be summed up only over a discrete set of values to test for 
stability, we can get a negative value of U2 with different 
simple elections. One way is choosing a linear function 
U {r) ^ a — (3r [fi > 0). In order to get a reasonable 
potential we have to cut it off at large distances. Beyond 
r = ri = a/P the potential would be taken as zero, and 
for r < tq (< ri) a strong hard-core will be supposed to 
avoid a complete collapse of the particles. This interac- 
tion will be referred to as the hard-core plus linear-ramp 
potential. B For this potential, U2 is negative as long as 
there is no particles at distances ro or ri from each other. 
In particular, if we restrict to values of riand ro such that 
fi/ro < we can conclude that a TS with lattice pa- 
rameter a, with ro < a < ri is unstable. The question 
is, what is the MEC in this case? A possibility is that 
at densities where the TL (taken as stable) has a lattice 
parameter a between ro and ri , the system segregates in 
two parts, both triangular lattices with lattice parame- 
ters ro and ri. In terms of the external pressiire P this 
would correspond to an isostructural transitiono at some 
pressure. If an isostructural transition exists, it means 
that the energy as a function of the volume of the system 
has a region with negative second derivative. Since we 
are considering changes of volume that do not change the 
crystalline structure, we can derive this necessary condi- 
tion for the existence of an isostructural transition easily 
from expression ([|), and the result is 
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< 0. 



(6) 



This condition is not satisfied by the hard-core plus 
linear-ramp potential. The conclusion is that in some 
range of pressures the TS is not the MEC of the sys- 
tem. The MEC for this potential, for different, values of 
P and ti/tq have been studied only recently^tl and will 
be presented in the next section. 



III. GROUND STATE FOR THE HARD-CORE 
PLUS SOFT REPULSIVE SHOULDER 
POTENTIAL 
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FIG. 2. Tfie family of potentials Ug (r). 

Here I will present the MECs for a family of potentials 
that vanish beyond some distance ri, are infinite for dis- 
tances lower that some tq , and in the intermediate range 
ro < r < ri are given by the expression 
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(7) 



which depends on the parameter g. When g ^ \ the 
potential reduces to a linear ramp between U {r^) — Uq 
and U {ri) = 0. For g ^ oo the potential has a square 
shoulder of height Uq between rg and ri , and for g — > it 
reduces to the simple hard-core potential at r = ro . The 
form of the potential for different values of g is shown 
in Fig. H This potential for particular values of g, and 
other related potentials have been studied since many 
years ago p^B but usually with the idea of the isostructural 
transition in mind, and part of the richness of the model 
has been missed (however, see also [7]). 
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FIG. 3. Minimum energy configurations for particles inter- 
acting through the potentials Ug (r), in the g-P* plane for 
ri/ro = 1.5 (a), and in the ri/ro-P* plane for g — > cxd (square 
shoulder potential) (b). Dotted lines would be the locations 
of isostructural transition between triangular lattices if other 
structures did not exist. 

The MECs for a fixed value of ri/rg = 1.65 as a func- 
tion of g and P* , and for a fixed value of g 00 (square 
shoulder potential) as a function of ri/ro and P* are 
shown in Fig. || (the adimensional pressure P* is defined 
as P* = rgP/Uo). The variety of configurations is in- 
triguing, and some of them could be guessed only after 
doing some numerical simulations, annealing from high 
temperature configurations. Note that a direct transi- 
tion between two TSs (one close packed, and the other 
expanded) is always preempted by the existence of addi- 
tional lower energy structures. The structures in Fig. ^ 
were found by inspection, and they are the lowest energy 
configurations found within each region, but other (more 
stable) structures may have been missed. Note that some 
of the structures have more than one atom per unit cell 
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(2 for 5*3 and 5 for S4) and there may be unequivalent 
sites within the structure (2 for 5*4). 

It can be mentioned here the interesting fact that for 
potentials with a positive value oi U2 only at a discrete 
set of values of r, the compressibility of the system (at 
zero temperature) is zero, or the system is anisotropic (in 
the sense that second order tensors are not necessarely 
proportional to the identity). This is due to the above 
mentioned fact that stability of the structure requires the 
existence of particles at distances at which U2 is positive. 
If the system is isotropic (symmetry C3 or higher) an in- 
finitesimal change in volume would make the number of 
particles at these distances be zero, and the structure 
destabilizes. So if the compressibility is different from 
zero, then the structure can have only a rotational sym- 
metry C2 (as it happens for instance with structure Sq 
in Fig. ^ , and the structure continuously deforms under 
changes of pressure, always keeping particles at distances 
where U2 is positive. 



IV. STABILITY AGAINST GLOBAL 
DEFORMATIONS: SCREENED CHARGES IN 
TWO DIMENSIONS 

The stability against displacements of a single particle 
is by no means sufficient to guarantee the global stability 
of the TS. Let us consider another kind of perturbation 
of the TS, consisting in a rescaling of all x coordinates 
of the particles by a factor p, and all y coordinates by 
a factor p~^, as illustrated in Fig. ^. This transforma- 
tion preserves the volume per particle in the system so 
stability is achieved if the energy has a minimum around 
p = 1. For p very close to 1 we can take p = 1 + e, 
with e <C 1, and do an expansion of the energy around 
e = up to second order. The result for the energyO is 
E = Eo + f/2£V2, with 



U2 



E 



U" in) + 3 



(8) 



may be 



For potentials with U' (r) < 0, expression 
negative even when (^ is positive. 

For instance, for inverse power interactions U (r) 
Uq (r/ro)~^ , expression (||) takes the form 



7 < 2, expression (g) diverges, and a correct calculation 
of U2 should take into account a long distance cut-off 
of the interaction (or the existence of the edges of the 
system). The function /2 (7) , which allows to calculate 
U2 using (^ , is shown in Fig. |l]. For 7 > 2 it is calculated 
directly from expression (llOh. For 7 < 2 it is calculated 
using (|^), with an exponential cut-off in U (r) of the form 
~ exp(— r/ro), with ro —^ 00. Note however, that the 
particular form of the cut-off does not influence the value 
obtained for U2- As we can see from the figure, for all 7 
the TS is stable. As for U2, U2 vanishes when 7^0. 




FIG. 4. Rescaling of tlie TS by a factor p along the x di- 
rection and in tlie y direction. 

It is interesting to note that for 7 < 2, the contribu- 
tion to U2 from any particle is negative, only the exis- 
tence of the cut-off makes the result to be positive. We 
can gain some insight on this point by a particular ex- 
ample. Let us take 7=1, and a sharp cut-off at some 
distance rg. To avoid the existence of sharp edges in 
the potential we shift the repulsive part so as to van- 
ish at ro, i.e, we will consider a potential of the form 
U {r) = 6' (ro — r) (l/r — l/ro), and calculate the energy 
as a function of p for different values of ro. We see (Fig. 
^ that although d^E/dp^ is always negative at p = 1, 
the interval around p = 1 with this characteristic narrows 
when ro increases, and in the limit of very large ro, the 
value of d^Ejdp^ becomes positive at p = 1 as soon as 
we smooth the cut-off. 



with 



U2 - h (7) f/o (r-o/a)' 



/2(7) 



7(7-2) 



(9) 



(10) 



which is only positive for 7 > 2. For 7 < 2 the negative 
sign of U2 would suggest that the TS is unstable. Note 
however, that 7 = 2 coincides with the value below which 
\J2 is dominated by long range interactions. In fact, for 
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FIG. 5. Energy per particle of the TS as a function of 
the rescaling parameter p, for particles interacting through 
a cut-off interaction, for different values of the cut-off ro. 
The curves were vertically displaced to make them coincide in 
p — 1. Note in the inset that for all ro the value of d^E/dp^ 
at p = 1 is negative, so the TS is unstable. However, the 
instability range rapidly diminishes with ro. 

The last example raises the follovifing question: what 
is the minimum sharpness of the cut-off necessary to get 
a negative value of U2 for some range of values of the 
cut-off parameter ro? For the family of potentials of the 
form exp [— (r/ro)'^] /r''', a numerical calculation based 
on (H) shows that for /i < 2, U2 is positive for any rg and 
7. For /J, > 2 and if 7 is lower than some value 70 (/i), 
there exists a range for ro such that U2 is negative. The 
function 70 (/i) goes to zero when /i — > 2+, and increases 
with fi. It becomes 1 for /i slightly larger than 3. In 
all cases the instability region occurs when the cut-off 
parameter ro is close to the lattice parameter of the TS. 



V. SUMMARY AND FURTHER DISCUSSIONS 

It was shown in this paper that identical particles inter- 
acting through repulsive central forces in two dimensions 
may have a minimum energy configuration very differ- 
ent than the usually expected triangular structure. This 
may happen for interaction potenticik as simple as the 
hard-core plus linear-ramp potential,EI or for sharply cut 
off power law interactions of the form 6 (ro — r) r ' , if 
7 < 2. 

We have concentrated here in the two-dimensional 
case. But it is not difficult to see that some of the conclu- 
sions can be extended to three dimensions. For instance, 
it is immediate to generalize expressiO|n_(H) to three di- 
mensions. The equivalent expression islij 
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Also the necessary condition ^ for the existence of an 
isostructural transition can be generalized to: 
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r? < 0. 



(12) 



Expression (|11|) may be negative even if (^^ is not satis- 
fied. Again this happens, for instance, for the hard-core 
plus linear ramp-potential. In three dimensions the a 
priori expected structures for simple spherical potential 
are the high symmetry structures sc, bcc, fee, or hep. 
Are they the only MFCs of this potential? The answer 
is negative. In Fig. |^ we see a preliminary diagram of 
MECs for this potential as a function of P* (now defined 
as P* = rgP/Ua) and ri/ro. This was obtained searching 
for the MFC among all crystalline systems with no more 
than two parameters determining their structure (this re- 
striction was used only to facilitate the search). These 
are the cubic (including sc, bcc, and fee Bravais lattices), 
tetragonal (simple (t) and body centered (bet)), rhom- 
bohedral (rh), and hexagonal systems. Only structures 
with one atom per unit cell were considered for simplic- 
ity, with the only exception of hexagonal structures (h) , 
where the closed packed structure (hep, with two atoms 
per unit cell) was included considering its well known 
stability. 
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FIG. 6. MECs for the hard-core plus linear-ramp potential 
in three dimensions. The search was performed among struc- 
tures of the cubic, tetragonal, rhombohedral, and hexagonal 
systems. Dashed regions are zones where none of these can 
be the MEC. See the text for more details. 

The results of Fig. ^ show that almost all of these 
structures are the MEC among the considered ones in 
some region of the ri/rg-P* plane. In addition, and hav- 
ing in mind the configurations found in the two dimen- 
sional case, it would not be surprising that other struc- 
tures corresponding to lower symmetry crystalline sys- 
tems, or structures with a more complex unit cells exist. 
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In connection with this, notice that the dashed regions in 
Fig. ^ correspond to lowest energy configurations (among 
the ones already mentioned) that have no particles at dis- 
tances ro or ri , and according to previous discussions we 
know that this structure cannot be stable, so the MEC 
is none of the considered ones. Clearly, numerical simu- 
lations are needed to exhaustively find all MECs for this 
or other related potentials. 

Another point that was not touched upon in this paper 
is the problem of stability at finite temperatures. For the 
hard-core plus linear-ramp potential in two dimensions a 
detailed discussion has been given elsewhere.B I will only 
mention here the interesting fact that in some cases the 
melting of the crystalline structures is anomalous (in the 
sense that it occurs with an increasing in density) due to 
the sudden availability of configuration space at higher 
energies upon melting (which may be assimilated to an 
effective reduction of particle size at melting). 

Other interesting issue concerns the dynamics of these 
structures. For instance, for the hard-core plus linear- 
ramp potential (at zero temperature), if we increase the 
external pressure smoothly, there is a value at which the 
TS would have a lattice parameter lower than ri, and 
the structure destabilizes against displacement of single 
particles. Since particles in different positions will move 
in rather independent directions, we expect to obtain a 
disordered (metastable) structure at high pressures. For 
potentials such as the sharply cut off 1/r, the instability 
is against deformations involving a large number of par- 
ticles, and thus the metastable structures obtained by in- 
creasing pressure are expected to consist of large patches 
of particles, deformed along different directions. This is 
in fact what is obtained in numerical simulations. This 
phenomenon, as well as the appearance of metastable 
structures when decreasing the temperature from a fi- 
nite value to zerou are important in connection with the 
transition to the glass state. 

Although the kind of potentials needed to obtain the 
behaviors discussed in this work are difficult to find in 
atomic systems, it is reasonable to expect that they have 
physical realizations in colloidal dispersions JiJ where the 
interaction potential between particles can be changed a 
great extent through the applications of different tech- 
niques. 
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